
 #For convenience we import the sample routine (so we can write sample(...) instead of
 #Bio::Phylo::EvolutionaryModels::sample(...).
 use Bio::Phylo::EvolutionaryModels qw (sample);
 
 ################################################################################
 #Simulate a single tree with ten species from the constant rate birth model with parameter 0.5
 ################################################################################
 my $tree = Bio::Phylo::EvolutionaryModels::constant_rate_birth(birth_rate => .5, tree_size => 10);
 
 ################################################################################
 #Sample 5 trees with ten species from the constant rate birth model using the b algorithm
 ################################################################################
 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              algorithm => 'b',
                              rate => 1,
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth,
                              model_options => {birth_rate=>.5});

                              
 #Print a newick string for the 4th sampled tree                              
 print $sample->[3]->to_newick."\n";            
 
 ################################################################################
 #Sample 5 trees with ten species from the constant rate birth and death model using 
 #the bd algorithm and two threads (useful for dual core processors)
 #NB: we must specify an nstar here, an appropriate choice will depend on the birth_rate
 #    and death_rate we are giving the model    
 ################################################################################
               
 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              algorithm => 'bd',
                              rate => 1,
                              nstar => 30,
                              threads => 2,
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth_death,
                              model_options => {birth_rate=>1,death_rate=>.8});
                               
 ################################################################################
 #Sample 5 trees with ten species from the constant rate birth and death model using 
 #incomplete taxon sampling
 #
 #sampling_probability is set so that the true tree has 10 species with 50% probability,
 #11 species with 30% probability and 12 species with 20% probability
 #
 #NB: we must specify an mstar here this will depend on the model parameters and the 
 #    incomplete taxon sampling parameters
 ################################################################################
                   
 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              algorithm => 'incomplete_sampling_bd',
                              rate => 1,
                              nstar => 30,
                              mstar => 12,
                              sampling_probability => [.5, .3, .2],
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth_death,
                              model_options => {birth_rate=>1,death_rate=>.8});

 ################################################################################
 #Sample 5 trees with ten species from a Yule model using the memoryless_b algorithm
 ################################################################################
 
 #First we define the random function for the shortest pendant edge for a Yule model
 my $random_pendant_function = sub { 
 	%options = @_;
 	return -log(rand)/$options{birth_rate}/$options{tree_size};
 };
 
 #Then we produce our sample
 my ($sample,$stats) = sample(sample_size =>5,
                              tree_size => 10,
                              algorithm => 'memoryless_b',
                              pendant_dist => $random_pendant_function,
                              model => \&Bio::Phylo::EvolutionaryModels::constant_rate_birth,
                              model_options => {birth_rate=>1});

 #################################################################################
 #Sample 5 trees with ten species from a constant birth death rate model using the 
 #constant_rate_bd algorithm
 ################################################################################

 my ($sample) = sample(sample_size => 5,
                       tree_size => 10,
                       algorithm => 'constant_rate_bd',
                       model_options => {birth_rate=>1,death_rate=>.8});

 #################################################################################
 #Sample 5 trees with ten species from a constant birth death rate model using the 
 #constant_rate_bd algorithm
 ################################################################################

 my ($sample) = sample(sample_size => 5,
                       tree_size => 10,
                       algorithm => 'constant_rate_bd',
                       model_options => {birth_rate=>1,death_rate=>1});

print "\n".'BD'.$sample->[0]->to_newick();                       
                       